Persistent supersolid phase of hard-core bosons on the triangular lattice 
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We study hard-core bosons with unfrustrated hopping (t) and nearest neighbour repulsion (U) 
on the triangular lattice. At half-filling, the system undergoes a zero temperature (T) quantum 
phase transition from a superfluid phase at small U to a supersolid at U c ~ 4.45 in units of 2t. This 
supersolid phase breaks the lattice translation symmetry in a characteristic y/% x \/3 pattern, and 
is remarkably stable — indeed, a smooth extrapolation of our results indicates that the supersolid 
phase persists for arbitrarily large U/t. 
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Introduction: The observation of strongly correlated 
Mott insulating states and T = superfluid-insulator 
transitions of ultracold bosonic atoms subjected to opti- 
cal lattice potentials^ has led to a great deal of interest in 
strongly correlated lattice systems that can be realized in 
such experiments^^ The recent observation of a super- 
solid phase in Helium 4 leads, in this context, to a natural 
question: Can the lattice analog of this, namely a super- 
fluid phase that simultaneously breaks lattice translation 
symmetry, be seen in atom-trap experiments? 

One class of promising candidates are systems which 
are superfluid when interactions are weak, but form insu- 
lators with spatial symmetry breaking when interactions 
are strong: In terms of conventional Landau theory, a 
direct transition between these two states is generically 
either first order, or pre-empted by an intermediate su- 
persolid phase with both order parameters nonzero; both 
types of behaviour are known to occur in specific lat- 
tice modelsi^Sii Moreover, as has been shown recently 
by Senthil et. al.fi conventional Landau theory can fail 
in certain situations in which quantum mechanical Berry 
phase effects produce a direct second-order phase transi- 
tion, thereby ruling out an intermediate supersolid phase. 
When such a transition occurS)2ii£ it is associated with 
quasi-particle fractionalization and deconfinement, and 
this alternative to an intermediate supersolid phase is 
thus interesting in its own right. 

Bosons on the triangular lattice with on-site repulsion 
V, repulsive nearest neighbour interaction U, and unfrus- 
trated hopping (t) provide a particularly interesting ex- 
ample in this context since the structure of interactions 
is simple enough that it can be realized in atom-trap 
experiments^iii In the hard-core V — > oo limit (which 
is also experimentally feasibl e - 2 ! 11 ) this maps to a sys- 
tem of S = 1/2 spins (Sf = rii — 1/2 where m is the 
boson number at site i) with antiferromagnetic exchange 
J z = U between the z components of neighbouring spins, 
ferromagnetic exchange J± = 2t between their x and y 
components, and magnetic field in the z direction equal 
to the chemical potential \i. It is this hard-core limit we 
consider in some detail below at zero chemical potential. 

Clearly, the ground-state in the limit U/t — > must 



be a featureless superfluid. On the other hand, the in- 
teraction energy U dominates in the U/t — > oo limit and 
leads to frustration since it is impossible to have all pairs 
of neighbouring spins pointing anti-parallel to each other 
along the z axis on the triangular lattice. The ground 
state in this limit is thus expected to live entirely in 
the highly degenerate minimally frustrated subspace of 
configurations with precisely one frustrated bond (par- 
allel spins) per triangle, and is selected by the dynam- 
ics associated with the hopping term t. The minimally 
frustrated subspace can be conveniently represented by 
noting that each state in this subspace corresponds to a 
perfect dimcr cover of the dual hexagonal lattice (with 
every frustrated bond on the triangular lattice mapping 
to a dimer placed on the dual link perpendicular to the 
bond in question). In this language, the effective Hamil- 
tonian in the U/t — > oo limit is then a quantum dimer 
model with a ring-exchange term that operates on each 
pair of adjacent hexagons of the dual lattice (see Fig.QJ. 

Quantum dimer models with ring-exchange on individ- 
ual plaquettes of two dimensional bipartite lattices quite 
generally have crystalline ground states in which the spa- 
tial arrangement of dimers breaks lattice symmetry>i2ii^ 
In our problem, a wavefunction that gains kinetic en- 
ergy from the double-hexagon ring exchange process on 
a maximal set of independently flippable hexagon pairs 
(see Fig.^l provides a similar candidate lattice symmetry 
breaking insulating state at large U/t (see Fig. 0. Note 
however that this analogy to simpler quantum dimcr 
models misses the important U(l) symmetry associated 
with charge conservation. Alternatively, one can focus 
on this U(l) symmetry at large U ft by thinking in terms 
of superfluid wavefunctions projected into the minimally 
frustrated subspace: Clearly, superfluidity can survive 
in such a projected state since the minimally frustrated 
manifold admits considerable charge fluctuations, and 
such wavefunctions also provide substantial kinetic en- 
ergy gain^ The breaking of lattice translation symme- 
try that seems natural by analogy to the simpler quan- 
tum dimer models then motivates a large- U variational 
ground state obtained by projecting a supersolid wave- 
functioni 1 ^ This suggests that the 'intermediate' super- 
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solid phase of Landau theory may, in fact, persist to large 
U in this case (another argument for a supersolid was 
given in Ref. ITfih . 

While these considerations are not definitive, they do 
at least emphasize that the behavior of this system at 
intermediate and large U presents very interesting pos- 
sibilities, and a detailed numerical study is one way to 
decide between them. In the present work, we use Quan- 
tum Monte-Carlo (QMC) methods to perform such a nu- 
merical study. Our results for the different T = phases 
are shown in Fig. ^ To summarize, we find that the 
superfluid at small U undergoes a transition to a super- 
solid phase at U c ~ 4.45 in units of 2i. This supersolid 
phase breaks lattice translation symmetry in a charac- 
teristic y/3 x v3 pattern shown in Fig. and appears to 
be indeed stable for arbitrarily large values of U / 1 (albeit 
with progressively smaller superfluid density). 

Model and method: Our Hamiltonian reads 



H = Mm ~ V2) (nj - 1/2) -t(b\ bj + hbt) 



(ij) 



+ J2[V(n i -l/2) 2 - f m i ] 



(1) 



where (ij) refer to nearest neighbour links of the two- 
dimensional triangular lattice, rii is the particle number 
at site i, and b\ is the boson creation operator at site i. In 
this work, we take the limit V — > oo to enforce the hard- 
core constraint, thereby mapping it to the S — 1/2 spin 
model as mentioned earlier, set t to 1/2, and take /i = 0. 
We use the well-documented stochastic series expansion 
(SSE) QMC method 16 i 17 i 18 which efficiently samples the 
high-temperature expansion of the grand-canonical parti- 
tion function. (At large values oiU/t, some modifications 
to the standard algorithm are necessary, and these will 
be discussed separateljiiS). 

Most of our data is on L x L samples with periodic 
boundary conditions and L a multiple of six ranging 
from 12 to 48 at inverse temperatures (3 ranging from 
10 to 30. Our choice of boundary conditions and as- 
pect ratio ensures that all the lattice symmetries are 
preserved after imposing the boundary conditions (see 
Fig. ^| . The nature of the T — ► phase and its low en- 
ergy spectrum of excited states is conveniently character- 
ized by the superfluid density p s , and the momentum (q) 
and imaginary time (r) dependent correlation functions 

Cp(q,r), C^f (q, r) of local particle density and kinetic 
energy respectively (a and a refer to the three possi- 
ble bond orientations T /i/ 2 shown in Fig. QJ. We use 
standard SSE estimators 18 to calculate p s , C p (q,T = 0), 
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J P drC p (q,r) 
r0 



C&%t = 0) and 



S p (q,Lj n 

S^(q,u> n — 0) = Jq drCK(q,T). These momentum 
space correlation functions are an unbiased probe of spa- 
tial order in the system. By analyzing the L and /3 de- 
pendence of the Bragg peaks at ±Q = ±(27r/3, 27r/3) 
seen in the equal time and static correlation functions of 
density and kinetic energy, we conclude that spatial order 
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FIG. 1: a) A flippable pair of spins and the corresponding 
hexagon-pair with its flippable dimer configuration. Also 
shown are static spins that surround this flippable pair to form 
an elongated hexagonal tile. Tiling the lattice with these tiles 
gives a candidate insulating state at large U/t (not observed 
numerically), b) Actual T — phase diagram and nature of 
spatial symmetry breaking in the supersolid phase. Darker 
bonds and sites indicate higher kinetic energy and density re- 
spectively, and the state shown corresponds to 6k — On = 0. 
Note that the lattice is drawn to emphasize periodicity in 
directions To and T2. 
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FIG. 2: Superfluid density p s and density wave order param- 
eter m 2 = S p (Q,uj n = 0)//3L 2 extrapolated to T — > and 
L — > 00. 



is established at these wavevectors when lattice transla- 
tion symmetry is broken in the supersolid phase (in the 
convention used above, the components of Q refer to pro- 
jections in directions To and T2 shown in Fig.^) 

We also measure two complex order parameters sensi- 
tive to this spatial symmetry breaking to obtain a better 
characterization of the supersolid state. These are de- 
fined as 



n A + n B e 2 "/ 3 + n c e 4 " /3 , 



1>k = K a + K h e 2 "/ 3 + K c e im / 3 



(2) 



Here the subscripts refer to the three-sublattice decom- 
position of the triangular lattice into A, B, C type sites, 
and a = BC, b = CA, c = AB type bonds respectively, 
while n and K are the densities and kinetic energies on 
the corresponding sites and bonds, ipx may be obtained 
from the Fourier components of the kinetic energies in 
the three lattice directions, \ using the relation 

IpK = e 4-/3^(°) + #*i/3 K M + e 4W3^(2) ; while ^ 
WWW 

is precisely equal to n^, the Fourier component of the 
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FIG. 3: Extrapolations implicit in Fig. |5] shown here for two 
values of U 
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FIG. 4: Binder cumulant g = 1 — (m 4 }/3(m 2 ) 2 as a function of 
U in the transition region. From the location of the crossing 
point we identify a T — phase transition to the supersolid 
phase at U « 4.45. Inset: Histogram of |i/>„| has a single peak, 
indicating a second-order transition 



density at the ordering wavevector Q. We expect both 
order parameters to average to zero as long as our al- 
gorithm remains ergodic — the probability distribution of 
their phases Ok and n however provides useful informa- 
tion regarding the nature of the supersolid phase. 

Numerical results: Our numerical results for the vari- 
ation in the superfluid density p s as a function of U are 
shown in Fig. El Each point shown in Fig. [21 represents 
an extrapolation of available data to the T = thermo- 
dynamic limit. The results summarized in Fig. [2] show 
no indication of any finite U/t quantum phase transition 
beyond which p s may become zero at zero temperature. 
Indeed, a smooth extrapolation of our data suggests that 
superfluidity persists in the low temperature limit at all 
finite values of U, albeit with an increasingly small T = 
value of p s . While this is surprising from the perspective 
of putative Valence Bond Solid (VBS) ground states of 
the corresponding quantum dimer model, further insight 
can be obtained by performing a variational calculation 
using projected superfluid wavefunctions; this work will 
be reported separately^ 

Although superfluidity survives in the entire range of 
U studied, the state at small U is not continuously con- 
nected to that at large U. Indeed, we see clear evidence 
for a continuous T — quantum phase transition at 
U c ~ 4.45. This transition point is estimated using stan- 
dard criteria in terms of Binder cumulants as shown in 
Fig. 0J (further details regarding the phase transition 
will be reported separate!^) . For U > U Cl the sys- 
tem spontaneously breaks lattice symmetry to produce a 
supersolid phase. To understand the nature of the super- 
solid phase, it is useful to analyze the joint probability 
distribution of the phases Ok and n defined earlier. At 
U = 10, we see from Fig. that Ok is essentially pinned 
to be equal to —20 n (modulo 2n) at low temperature 
and large L, while 9 n has a distribution which peaks at 



9f t = 2irp/6 with p an integer from to 5. The pic- 
ture that emerges (and gets sharper at larger L and /?) 
is thus of a state in which a relatively more mobile fluid 
of density phx living on a hexagonal lattice backbone of 
the full lattice is responsible for the superfluidity, while 
the centers of the hexagons have an average density of 
p c that is less mobile. The six values of p correspond to 
three possible hexagonal backbones of a triangular lattice 
in conjunction with two choices for the sign of phx — Pc- 
Note that this spatial order is accompanied by a very 
slight deviation of the total density p from 1/2 (which 
survives in the T = thermodynamic limit), with the 
sign of deviation given by that of phx — Pc- We have also 
monitored these histograms at larger U < 40. We find 
that Ok remains strongly pinned to —28 n , and while the 
pinning of 9 n and 9k individually does become weaker, 
the basic picture of the supersolid state remains the same. 

Landau theory: Much of this picture of the supersolid 
phase can be understood within the framework of a Lan- 
dau theory written in terms of the order parameters ip n 
and tpK (while it is not necessary to do so, we find it 
convenient to explicitly include ipK in our description). 
For our purposes here, it suffices to consider only the 
'potential energy' terms of the Landau theory and leave 
out all fluctuation terms that involve spatial and time 
derivatives, or couplings to the superfluid order parame- 
ter, although these can also be straightforwardly written 
down. Terms in the Landau theory are constrained by 
the requirement of invariance under all the symmetries 
of the system. The action of these on our order parame- 
ters is simple to state: Under both lattice translations To 
and T 2 we have ^ n -> e 2 "/ 3 ^„, -> e 27ri ^ifj K , while 
under a 7r/3 rotation about a A sublattice site, we have 
tpn —> ipni ^ K ~ * ^P*k ■ Finally, ip n is odd under particle- 
hole transformations, while ipK is even. Terms consistent 
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(a) U = 10, 13 = 30, L= 18 o 501 



(b) — U=10 |S=20 
U = 40 11 = 20 




FIG. 5: Top panels: Greyscale plot of the joint probability 
distribution of 9 n and 9 k', 9 n dependence of p. Bottom pan- 
els: Probability distribution of 9„ 9 k + 20 n , and 9k- The his- 
togram of 9k + 20 n has additional peaks at ±27r (not shown). 



with these symmetries at pi = give, up to sixth order 

S po t&n, 4>k) = /(|Vn| 2 , \M 2 ) + CoM + ^n) 

+C0k(^K + i>*K ) + CnK^K + C^k) ■ (3) 

As usual, spatial symmetry breaking corresponds to the 
function S po t developing a minimum at a nonzero value 
of \4>n\- The detailed nature of the ordering is determined 
by the signs of coefficients c which fix the relative as well 
as absolute phases of the two order parameters. The 
results shown for U = 10 in Fig. [5] can be modeled by 
taking all c negative, and this translates to the schematic 
picture of the phase shown in Fig.^ In addition, the very 
slight 9 n dependence of p— 1/2 can be modeled^ by the 
presence of a coupling term (p — l/2)(^ + V^ 3 ) with a 



tiny positive coefficient. 

Discussion: We have thus established the presence of 
a persistent low-temperature supersolid phase on the tri- 
angular lattice. The remarkable stability of this phase is 
in contrast to the relatively small window of parameters 
within which supersolids have been previously seen on the 
square latticeA^ Indeed, a smooth extrapolation of our 
data indicates that the supersolid phase persists in the 
U/t — > oo limit. Our results thus throw up an interesting 
possibility for the observation of this phase in atom-trap 
experiments. It would therefore be useful to map out the 
finite temperature diagram for U > U c . In the absence of 
any coupling between spatial and superfiuid order param- 
eters, superfluidity would be lost by a Kosterlitz Thouless 
(KT) phase transition, while crystalline order would be 
lost via two KT phase transitions with an intermediate 
power-law ordered crystal phased analogous to that seen 
in the transverse field Ising antiferromagnet on the same 
latticeiSiiSS The coupling between the two order param- 
eters is expected to modify the detailed nature of the 
finite temperature phase diagram; this will be reported 
on separately^ Our picture of the supersolid state in- 
dicates that it is not destabilized by doping, and this is 
currently under investigation. 19 Finally, it would also be 
interesting to study the stability of the supersolid upon 
relaxation of the hard-core constraint. 
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